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Abstract 



We discuss the direct use of cubic-matrix splines to obtain continuous approxima- 
jr \ • tions to the unique solution of matrix models of the type Y"{x) = f(x, Y(x)). 

For numerical illustration, an estimation of the approximation error, an algorithm 
for its implementation, and an example are given. 



1 Introduction. 

Matrix initial value problems of the form: 

;U : Y"{x) = f(x,Y(x)) 

QV } a<x<b, (1) 

Y(a) = Y , Y'(a) = Y t 

are frequently encountered in different fields of physics and engineering (see e.g. BZha02l ). 
In the scalar case, numerical methods for the calculation of approximate solutions of 
(•*) \ (Q]) can be found in [Col93|. For matrix problems, linear multi-step matrix methods 

with constant steps have been studied in HJod93l . Although in this case there exist 
a priori error bounds for these methods (expressed as function of the data problem), 
these error bounds are given in terms of an exponential which depends on the integra- 
tion step h. Therefore, in practice, h will take too small values. Problems of the type 
(fl} can be written as an extended first-order matrix problem. Such a standard approach, 
however, involves an increase of the computational cost caused by the increase of the 
problem dimension. Recently, cubic-matrix splines were used in the resolution of first- 
order matrix differential systems [Def06|, obtaining approximations that, among other 
advantages, were of class C 1 in the interval [a, b], and easy to compute producing an 
approximation error 0(h A ). The present work extends this powerful scheme to the so- 
lution of matrix problems of type ([T). Throughout this work, we will adopt the notation 
for norms and matrix cubic splines as in |Def06| and common in matrix calculus. The 
paper is organized as follows. Section [2] develops the proposed method. Finally, in 
Section[3j an example is presented. 
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2 Construction of the method. 

Let us consider the initial value problem 

Y"(x) = f(x,Y(x)) 

a<x<b, (2) 

Y(a) = Y Q , Y'(a) = Y x 

where Y 0l Y u Y{t) e C rxq , f : [a,b] x C rxq x i — >C rxq ,f € C° (T), with 

T = {(a>,Y); a<x<6, Y eC"«} , (3) 

and / fulfills the global Lipschitz's condition 

ll/(z,Fi) - /(x,y 2 )|| <i||yi-y 2 || ,o<x<6,Fi,y2ec- x «. (4) 

Let us also use the partition of the interval [a, b] defined by 

A[ a ,b] = {a = x < xi < . . . < x n = b} , x k = a + kh , k = 0, 1, . . . , n , (5) 

where h = (b — a)/n, n being a positive integer. We will construct in each subinterval 
[a + kh, a+(k + l)h] a matrix-cubic spline approximating the solution of problem©. 
For the first interval [a, a + h], we consider that the matrix-cubic spline is given by 

S \ la , a+h] W = *» + Y '{a){x-a) + ^Y"{a)(x-a) 2 + ^Aoix-a) 3 , (6) 

where A E C rX9 is a matrix parameter to be determined. It is straightforward to 
check: 

Si (a) = Y(a) , Si (a) = Y'(a) , S[' (a) - Y"(a) = /(a, S\ (a)) 

|[a,a+h] y ' y > ' \[a,a+h] V ' V ' ' | [a,»+h] W V ' V I [a,a+fc] V " 

Thus, © satisfies the equations of problem (f2j) at point x = a. To fully construct the 
matrix-cubic spline, we must still determine Aq. By imposing that © is a solution of 
problem §2& in x = a + h, we have: 

Si' (a + h) = f(a + h,S { +n (a + h)) , (7) 

and obtain from (0 the matrix equation with only one unknown matrix A$ : 

(8) 



*-i 



f(a + h, Y(a) + Y'{a)h + ^Y"(a)h 2 + ^A hA - Y"(a) 



Assuming that the matrix equation (O has only one solution Aq, the matrix-cubic spline 
is totally determined in the interval [a, a + h]. Now, in the next interval [a + h,a + 2h], 
the matrix-cubic spline is defined by: 

S\ (x) = S\ (a + h) + Si (a+h)(x-(a + h)) 

\[a + h,a + 2h] y > \[a,a + h] V ' |[a,o+h] V ' V V " 
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+ nS\ (a+h)(x-(a+h)f+-A 1 (x-(a+h)f , (9) 



so that S(x) is of class C 2 ([a,a + h] U [a + h,a + 2h)), and all of the coefficients of 
matrix-cubic spline S\ (x) are determined with the exception of A\ £ C rxq . 

By construction, matrix-cubic spline (O satisfies the differential equation (0 in x = 
a + h. We can obtain A\ by requiring that the differential equation (f2j) holds at point 

x = a + 2h: 



Si' 



(a + 2h) = f (a + 2h, S\ (a + 2h) 



\[a + h,a + 2h] 

Expanding, we obtain the matrix equation with only one unknown matrix A\\ 



f(a+2h,S\ (a+h)+S\ (a+h)h+-S? (a+h)h 2 + -A 1 h 3 ) 

{a+h)] . (10) 



- s \ 

\[a,a+h] 

Let us assume that the matrix equation ( TTOb has only one solution A\, This way the 
spline is now totally determined in the interval [a + h, a + 2h]. Iterating this process, 
let us construct the matrix-cubic spline taking [a + (k — l)h, a + kh] as the last subin- 
terval. For the next subinterval [a + kh, a + (k + l)h], we define the corresponding 
matrix-cubic spline as 

S \ [a+kh , aWW M = &(*) + k. Mx ~ {a + kh)f 

2 1 

where f3 k (x) = V -S, (0 (a + kh)(x - (a + kh)) 1 . (11) 

£ ' /! [a + (fc-l)h,a + fch] 

i=0 



With this definition, it is S(x) £ C 2 [ I) [a + jft, a + (j + l)/i] which fulfills the 

_ V=° / 

differential equation (f2]i at point a; = a + fc/i. As an additional requirement, we assume 

that S(x) satisfies the differential equation (O at the point x = a + (k + l)h, i.e. 
St' (a + (k + l)h) = f(a+(k+l)h,S\ (a+(k+l)h) 

\[a + kh,a + (k + l)h] V V ' ' \ ^ ' \[a + kh,a+(k + l)h] V V ' > 

Subsequent expansion of this equation with the unknown matrix A k yields 



*-s 



/ ( a + (k + l)h, fa{a +(k + l)h) + ^A k h 3 ) - ;,;'(„ - (A- - I )h i 



(12) 

Note that this matrix equation (fT2l is analogous to equations ([8]l and ( fTOb , when k = 
and k = 1, respectively. For a fixed /j, we will consider the matrix function of matrix 
variable g : C rxq ^ C rxq defined by 



9(T) 



f[a+(k + l)h, [3 k (a + {k + l)h) + -Th 3 ) - [3' k \a + {k + l)h) 



Relation ( fT2b holds if and only if Ak — g(Ak), that is, if Ak is a fixed point for function 
g(T). Applying the global Lipschitz's conditions (01, it follows that 

Lb? 
^(TxJ-flfCTa)!!^— ||Ti-T 2 || . 

Taking h < \/f-, g(T) yields a contractive matrix function, which guarantees that 
equation (fT2l) has unique solutions Ak for k = 0, 1, . . . , n — 1. Hence, the matrix- 
cubic spline is now fully determined. Taking into account 0Los67l Theorem 5], the 
following result has been established: 



Theorem 2.1 Ifh< \/ j-, then the matrix-cubic spline S(x) exists in each subinterval 

[a + kh, a + (k + l)h], k — 0,1, ... ,n — 1, as defined by the previous construction. 
Furthermore, if f € C l (T), then \\Y(x) — S(x)\\ = 0(h 3 ) \fx €E [a, b], where Y(x) is 
the theoretical solution of system (0. 

Depending on the function /, matrix equations (H) and (fT~2b can be solved explicitly or 
by using some iterative method HOrt72l . Summarizing , we have the following algo- 
rithm: 

• Take n > -j= , h = (b — a)/n and Au w defined by ©. 

• Solve dSl and determine Si (x) defined by ([6b. 

— \[a,a + h] X ' J 

• For k = 1 to n — 1, solve ( fT~2T >. Determine S\ (x) defined by (fTTT i. 

\[a + kh,a + (k + l)h] X ' J 

3 Example 

The problem 

Y"{t) + AY(t) = , (13) 

with Y(0) — Y , Y'(0) = Y t , has the exact solution 

Y{t) = cos (VAt)Y + (Va\ sin (VAt)Yi , 

where \f~A denotes any square root of a non-singular matrix A, IIHar05l . The prin- 
cipal drawback of this formal solution is the difficult computation of vG4, cos (yAt) 
and sin (y/At). The proposed method avoids this drawback. We consider problem 

O where A = ( \ \ ] ,Y - ( ° Q J ) ,^i = ( j J ), t € [0,1], whose 

1 



exact solution is Y(t) = sin 
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/-W .0 Yin this case 

tcos(i) sm(t) / 

L fa 2.82843. By Theorem I2TT1 we need to take h < 1.45647, so we choose ft = 0.1 
for example. The results are summarized in the following table, where the numerical 
estimates have been rounded to the fourth relevant digit. In each subinterval, we eval- 
uated the difference between the estimates of our numerical approach and the exact 
solution. The maximum of these errors are indicated in the third column. 



Interval 


Approximation 


Max. 


Error 


[0,0.1] 


/ x - 0.1664a; 3 
^ x - 0.4986k 3 


\ 

s - 0.1664s 3 / 






1.0072 


x 10~ 6 


[0.1,0.2] 


/ 1.00005x - 0.0005s 2 - 0.1647s 3 
^ 1.0002s - 0.0025s 2 - 0.4903s 3 




1.0001s - 0.0005s 2 - 0.1647s 3 


) 


6.3032 


x 10~ 6 


[0.2,0.3] 


/ 1.0005s - 0.0025s 2 - 0.1614s 3 
\ -0.0001 + 1.0022s - 0.0124s 2 - 0.4738 


\ 

s 3 1.0005s - 0.0025s 2 - 0.1614s 3 J 


2.0059 


x 10 -5 


[0.3, 0.4] 


/ -0.0002 + 1.0018s - 0.0069s 2 - 0.1565s 5 
^ -0.0008 + 1.0088s - 0.0344s 2 - 0.4494s 3 


-0.0002 + 1.0018s 



- 0.0069s 2 - 


0.1565s 3 J 


4.6213 


x 10 -5 


[0.4,0.5] 


/ -0.0006 + 1.0049s - 0.0147s 2 - 0.1500s 3 
^ -0.0028 + 1.0242s - 0.0728s 2 - 0.4174s 3 


-0.0006 + 1.0049s 



- 0.0147s 2 - 


0.1500s 3 J 


8.8359 


x 10 -5 


[0.5,0.6] 


/ -0.0016 + 1.0109s - 0.0266s 2 - 0.1420s 3 
^ -0.0077 + 1.0536s - 0.1316s 2 - 0.3782s 3 


-0.0016 + 1.0109s 



- 0.0266s 2 - 


0.1420s 3 J 


1.4964 


x 10~ 4 


[0.6, 0.7] 


/ -0.0036 + 1.0210s - 0.0436s 2 - 0.1327s 3 
^ -0.0176 + 1.1030s - 0.2140s 2 - 0.3324s 3 


-0.0036 + 1.0210s 




- 0.0436s 2 - 


0.1327s 3 J 


2.3267 


x 10" 4 


[0.7,0.8] 


/ -0.0073 + 1.0368s - 0.0661s 2 - 0.1219s 3 
^ -0.0354 + 1.1791s - 0.3227s 2 - 0.2807s 3 


-0.0073 + 1.0368s 




- 0.0661s 2 - 


0.1219s 3 J 


3.3941 


x lO" 4 


[0.8,0.9] 


/ -0.0134 + 1.0597s - 0.0947s 2 - 0.1100s 3 
^ -0.0646 + 1.2885s - 0.4595s 2 - 0.2237s 3 


-0.0134 + 1.0597s 




- 0.0947s 2 - 


0.1100s 3 J 


4.7114 


x lO" 4 


[0.9, 1] 


/ -0.0229 + 1.0914s - 0.1299s 2 - 0.0970s 5 
^ -0.1093 + 1.4378s - 0.6253s 2 - 0.1623s 3 


-0.0229 + 1.0914s 



- 0.1299s 2 - 


0.0970s 3 J 


6.2838 


x 10- 4 
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